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ABSTRACT 



WOLfYKHoxuBRAfir 



A new method of numerical solution for time-harmonic 
field scattering problems involving inhomogeneous bodies of 
revolution is presented. A simple single-loop multiport 
feedback system is driven by the incident field. The forward 
path transfer matrix, found by use of the finite-element 
method, translates the total near-field into the body-surface 
field. The feedback path matrix is established by the magne- 
tic vector potential formulation; it translates the body- 
surface currents into the scattered near-field. The thin- 
wire scattering problem is formulated in terms of this field 
feedback method,. Results of this computation are compared in 
graphical form to those found by Kallen's integral equation 
solution of the scattering problem for several cases of 
interest . 
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I. INTRODUCTION 



The objective of this endeavor is to investigate the per- 
formance of a new, and potentially quite powerful, technique 
for the computation of electromagnetic scattering problems. 

For reasons that will later become apparent, this new technique 

3 

will be denoted the "field feedback formulation" (FFF or F ) . 

3 

The work reported here is a first application of F , where 
the calculation of scattering from straight thin-wires has 
been performed. This is an almost classical electromagnetic 
field problem that has been routinely solved by use of inte- 
gral equations for many years. Such a cononical problem was 

3 

chosen as a testing vehicle for the newly conceived F 
analysis method since alternate integral equation solutions 
are readily available for comparisons. 

The formulation developed here employs a combination of 
the finite element numerical method and the classical magnetic 
vector potential. The finite element method (FEM) is an 
application of variational techniques which implements a 
weighted residual approximation. The high accuracy of the 
finite element method is used to advantage in reducing the 
amount of computing required, as compared to other numerical 
methods. Use is also made of coupled azimuthal potentials 
(CAP). The CAP formulation restricts the class of approachable 
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problems to those which are axially symmetric, but does allow 
the treatment of material inhomogeneities and multiwavelength 
dimensions . 

Thin-wire assumptions are made which greatly simplify 
the current-scattering problem. These assumptions, however, 
are by no means critical to the solution method. The formu- 
lation should be readily extendable to a large class of impor- 
tant interior, radiation, and scattering problems. 
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II. THE PROBLEM AND METHOD OF SOLUTION 



Electromagnetic scattering problems usually involve 
a known incident field, such as a uniform plane wave. The 
scattered field can be found from the current induced by the 
presence of the incident field, but the induced current can 
be found only after the total field is known. Since the 
total field is just the sum of the incident and scattered 
fields, a definite difficulty has been encountered. Know- 
ledge of the induced current is needed to find the scattered 
field, while at the same time, the scattered field must be 
known to determine the induced current . 

This is clearly seen in the form of Maxwell's equations 

for the scattering problem. Assume that a scatterer with 

properties Y^ = a^ + j ( jJS^ , Z^ = ja 3 y, is in free space where 

Y =joJU , Z =jwy , and in the presence of sources J , M . 
o-’ o o-' o ^ o’o 

Then the curl equations have the form 



VxH = J +Y,E 



o 1 



(la) 



VxE = -M -A,H 



(lb) 



o 



where the induced electric and magnetic currents are 



Jl -- (Y^-Y^)E 



(2a) 



M, = (Z,-Z )H 



1 o 
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and the scattered magnetic field is 



H 



scat 



= - H. 



incident 



( 3 ) 



The electric field relationship is, of course, similar. 

The simple and compact nature of the simultaneous equa- 
tions (1) and (2) is highly appealing, but of little practi- 
cal use in obtaining numerical answers to a problem because 
of the excessive computing required. Yet, while the require- 
ment to have two unknowns cannot be escaped, the equations 
relating these unknowns can at least be rendered more 
tractable. 

A. THE FIELD FEEDBACK FORMULATION 

As will be shown later, two simple matrix equations can 
be derived which relate the magnetic and electric fields at 
a number of points along the surface of the scatterer to the 
fields at points along a boundary some distance away. These 
equations are independent and work in opposite directions. 
That is, from the first equation, given the total boundary 
field values, e^ and h ^ , the scatterer surface total field 
values, e and h , can be found: 




(4a) 



T 

The column vectors, [e.'h.] , contain the field values at' the 

— 1—1 ’ 



I 



» T 

computation points. The second equation takes Ce I h ] and 



T 

yields the scattered field values [h ] along the boundary 

s I s 
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(4b) 



Ce I 



s I 



h ] 

s 



= T [e, Ih, ] 



— W —CO 



Equations (4a), (4b), and (3) comprise a simple feedback 
system, as in Figure 1. 




e 

— s 





Figure 1. Field Feedback 

What remains , then, is to construct the matrices V and T 
to fit the problem at hand. It should be pointed out that 
the magnetic field alone is involved in the wire-scatterer 
problem as developed here. So [e ih ] is replaced by h . 

The field translation matrices V and T, though developed in 
exactly the same way, will be considerably simplified as 
compared to the general scatterer case. 

B. THE WIRE SCATTERING PROBLEM 

Given a perfectly conducting thin wire in free space, and 
a plane wave of known amplitude incident from a known direction, 
it is required to find the current induced on the wire and the 
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field scattered by the wire at any point in space. Finding 
the induced current is equivalent to finding the field on the 

/s A 

wire surface, since J = n x H , where n is the unit vector 
normal to the surface. Normalized coordinates will be used 
throughout the problem, and are given by 

(R,Z,4>) = (K r, K z, <j>). 
o o 

The coordinates (r,Z,<j>) are standard circular cylindri- 
cal coordinates, and K =2ir/A is the free space wavenumber of 

o o 

the time-harmonic field. The situation is depicted in 

Figure 2. Use of these normalized coordinates will allow 

the formation of a compact and useful form of the CAP equations. 

The wire radius is a, given in terms of 1^, the free space 
wavelength. Then the normalized radius is R =27ra/X . The 
wire length is L, also given in wavelengths. So Z^=2'nL/X^. 

The incidence angle of the plane wave is a, measured as 
shown in Figure 2. Since perfect conductivity is assumed, 
fields and currents will exist only on the surface of the 
wire. And since a very thin wire will be assumed, where a 
<A^/10, the current at the wire ends may be set to zero. 

The problem will not be solved in terms of the usual 
electric and magnetic fields, E and H, or the scalar poten- 
tial, V. Instead, use will be made of the magnetic vector 
potential, A. The current and scattered fields will be 
presented and calculated in terms of the CAP formulation 
discussed below. 
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Figure 2. Normalized Coordinates 

C. COUPLED AZIMUTHAL POTENTIALS 

The CAP formulation [3] has been developed to handle 
previously intractable time-harmonic field problems involving 
continuously and discretely inhomogeneous bodies of revolu- 
tion. With this limitation of axial symmetry, all electro- 
magnetic field components can be expanded via exponential 
Fourier series in the <p coordinate. The model fields can 
then be represented by two coupled azimuthal potentials . 

These potentials satisfy a system of coupled partial differ- 
ential equations, as well as a variational criterion, which 
will be developed here. 

Using the normalized coordinates in Figure 2, the total 
fields are decomposed into azimuthal modes: 
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E(R,Z,(f>) 



(5a) 



00 




(R ,Z ) exp ( jm(}) ) 



n^HCRjZ ,<()) 



00 

( R , Z ) exp ( j m(j) ) 

m=-oo 




(5b) 



where n^=120 7r. By substitution of these expansions into 
Maxwell’s equations, it can be shown that all modal field 
components can be generated from two coupled azimuthal 
potentials tj;^(R,Z,m) and i|; 2 (R 5 Z,m), using 

✓s /\ 

<J) X e m = jfj..^(m(|)xVii;^-Ru^Vii; 2 ) (6a) 



■ 



e 

m 






4) X h^ = jfj,^(m(j)xV4;2 + Re^V4;^) 



(6b) 
( 6c) 
(6d) 



where 4> is the unit azimuthal vector, V is the two dimen- 
sional gradient operator, and 

fj^[R,Z] = [y^(R,Z)e^(R,Z)R^-m^]'^ (7) 

The Euler-Lagrange variational technique will be dis- 
cussed more fully in the exposition of the finite element 
method below. It is sufficient to note here that this tech- 
nique finds the solution of the system partial differential 
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equations as the stationary point of a functional of the 
potentials and their first derivatives. This stationary 
point is defined as the point which the 

variation of F, 5F, vanishes for perturbations 6^^) 

about that point. 

The functional F is a surface integral of the Lagrangian 
L, where 

The functional F is integrated over the constant azimuth 
planar meridian cross section Q of the volume in which the solu- 
tion is of interest. 

It is convenient at this point to change variables , and 
introduce a simplifying assumption. From (6b) and (6d) we 
have e^^(R , Z ) =i|;^/R and h^^(R,Z) = So finding only the 

azimuthal components of the modal fields, e^^ and will 

yield a complete solution through (6a) and (6c). Now, since 
a thin wire is assumed, the phase variation of the incident 
field from side to side on the wire may be neglected. Then 

/s 

the induced current will be a function of Z alone; 

Therefore the scattered fields will be azimuthally invariant, 
so that E(R,Z,(1)) = e^(R,Z) and p^H(R,Z ,cl))=h^(R,Z) . Further- 

A 

more, since J„=nxH=H,Z and since 
Z cp 
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1 

j we 



1 

j we 



VxH 



A 

[-R 



. 2 1 A 

az R SR 



R H^], 



it follows that e,^ = 0, and h = h, . The number of modes 

(po o <po 

has now been reduced to one, the h, mode. 

(J)0 

Observe that (6a) and (6c) no longer involve coupled 
potentials. This allows the Lagrangian to be decoupled, and 
it may now be written 



= f {V(Rh.^) 

O (j)0 



[R ■V(Rh^ )]}-y Rh: 
yr (j )0 ^r ij)0 



So the functional is 



( 9 ) 



F 



/ 



L(R,Z,h, ,Vh. )dRdZ 
’ (po (po 



(10) 



It remains to find h^^(R,Z) to make F stationary in for 
given Dirichlet boundary conditions . This will be done by 
use of the finite element method. 
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III. COMPUTATIONAL PROCEDURES 



The constructions of the V and T matrices have little 
in common other than the coordinate system. The finite 
element method is used to generate V. T is derived from 
the magnetic vector potential in this perfect conductor case, 
while if a dielectric body were used, both electric and 
magnetic vector potentials would be required to construct T. 

The use of numerical methods obviously implies that the 
current and fields will be found only at discrete points. 

As discussed in Section II. C., only points in a constant 
azimuth meridian plane need be considered. The region between 
the wire surface and the boundary at which the fields are to 
be calculated is divided into triangular elements. These are 
areas of integration in the finite element method. The divi- 
sions and the nodal numbering scheme are shown in Figure 3. 

There are a total of points along the wire, and Ng points 
along the boundary. 

A. FINITE ELEMENT METHOD 

This method is one whereby the region of interest is divided 
into a number of overlapping subregions of a specified shape. 
Then the problem solution function 9 is approximated by 
defining the solution as a linear combination of the products 
of unknown coefficients with "basis" functions defined over 
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over each of the overlapping subregions. These basis 
functions are usually sectionally linear, and, in very simple 
terms, are used to "fill in" the approximation between points. 
A weighted residual approximation of the solution 9, commonly 
known as the Ritz method, can then be made by the application 
of the discretized Euler-Lagrange variational procedure, 
subject to the appropriate boundary conditions. The rigorous 
development of each of these parts of the finite element 
method is not necessary here, but the background required for 
a basic understanding of the method will be given below. 



Z 








0 



Figure 3. Solution Points 
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1 . Basis Function Expansions 



In order to solve a problem numerically, some method 
of approximating the solution using a finite number of 
unknowns must be found. Also, the method’s convenience must 
be balanced against the desired accuracy, A general term 
for such an approximation is a "basis function expansion" 
of a function 0 [2]. The type of expansion used in the 



finite element method for a one-dimensional problem is given 
by 



where 6^^ is the discrete version of 6, the 6^ are the 
values of 6(X) at the discrete points X^, and the U^(X) 



functions are usually employed in FEM calculations, though 
quadratic or higher-order polynomial functions could be 
used to obtain greater accuracy. Linear basis functions for 
a one-dimensional problem are shown in Figure 4. Note that 
the node spacing is irregular. This allows tailoring of the 
expansion to better approximate the solution, and, in two 
or higher dimensional cases, may be necessary because of the 
shape of the solution region. If neither of these are 
necessary, regular node spacing can be used with some advan- 
tage to reduce complexity. The basis fucntions u^^ are seen 



N 




( 11 ) 



i = l 



are the basis functions at the points X^. Linear basis 
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to have unit value at their respective nodes , and are 
zero for ]_ 3 l_fi_fN . 




Figure 4. Linear Basis Functions 

A two-dimensional linear basis function often used is 
the "pyramid" function; its two-dimensional projection is 
shown in Figure 5. A rectangular region in the x-y plane 
has been discretized by division into the triangular elements 
shown. The boundary 3^22 ^ function U2 

been outlined to make the shape clear. 

The function has a value of one at the center node, and 
is zero along and outside the periphery. Its equation can 
be written 

U22(x3Y) = U32 j = ^12) 

where the component functions are planar segments of the form 

u ^2 - ax-B y+Y (13) 
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Figure 5. Two-Dimensional Finite Element Mesh 

Then a,S,Y must be functions of (x,y) so that u^ ^ *^^3 »y 2 ^ ’ 

u ^2 ^^i 5 Yj ) “0 » i?^3, j?^2, for integer (i,j). 

As in the one-dimensional case, the basis functions over- 
lap; that is, a function is centered around each node in the 
region or on the boundary. Of course, on the boundary of 
at least half of the function is "lost" outside Q. The 
overall appearance of the basis function array in three di- 
mensions might be imagined as a crystalline surface of multi- 
faceted peaks centered over the nodes, with raised valleys 
running between. If a function 9(x,y) were defined over , 
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its expansion could be written 



7 6 

i=3 i=i 

If the exact values of 0.. were known, 0,, would be a multi- 
faceted approximation to 0^^ , equal only at the nodes . 

Then, within the limitations of machine accuracy and computing 
time, the greater the node density the better the approxima- 
tion 0j^ will be . 

The wire scattering problem expands the field, h^^, by 

use of these pyramid basis functions. Referring to Figure 3, 
it is seen that the region covered by the triangular ele- 
ments will not contain any complete pyramids, but this is not 
really necessary. It will prove convenient to use a sequen- 
tial node coordinate system rather than the (i,j) Cartesian 
description of (12). Let n be the node coordinate, where 
l_<n_<N , N=N^+Ng. Then the approximating expansion h^^^ of h^^ 
can be written 

N 

n=l 



Note that the top and bottom wire nodes are not represented. 
This is because of the assumption that h^^=0 at these points. 

Two field vectors can now be defined. The nodal wire- 
surface field is 

h = h , n = 1, N (16) 

— w — n ’ ’ w 
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The nodal boundary field is 



= ^n ’ ^ N 



( 17 ) 



These vectors contain the approximate values of h^^ at the 
nodes of the triangular mesh. 

2 . Euler-Lagrange Variational Technique 

Consider the differential equation A9=f, which together 
with certain boundary conditions, describes some physical pro- 
blem in a region Q. The complex scalar wave (Helmholz) 

2 2 

equation, A =[V +k ]0=O is an example. The differential equa- 
tion is related to the functional F(9), where 



F(9) 




L(0, D 9 



D 9, 

y 



X 



y ) dQ 



(18) 



This relationship is a fundamental result of the calculus of 
variations, and simply stated, is the fact that the functional 
F(0) is stationary at the "point" 9^ which also satisfies the 
differential equation. That is, if the particular function 9^ 
can be found which causes the derivative (or first variation) 
of F(9) to vanish, then F(9) is said to be stationary, and 
0^ is said to be the stationary point of F(0); 0^ is then 
the solution of the equation A9=f [5]. 

Observe that the vanishing of dF/d9 implies that a 
minimum, maxim.um, or saddle point of F occurs 6 = 9^' In many 
physical problems F is associated with the potential energy 
of systems and thus finding 9^ is sometimes termed the 
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"minimization" of F. The problem now presents a choice of 
computational methods . The differential equation might be 
approximated by a discrete system, say finite differences, 
or the variational integral, F, can be "minimized" over a 
finite number of trial functions, as is done by the finite 
element method. 

It is necessary in this procedure to find the Lagrangian 
L for use in F(9) [4]. This L must be such that its insertion 
in the Euler-Lagrange equation produces the differential 
equation that is to be solved. As an example, the Euler- 
Lagrange equation for the two independent variable case is 
given by 
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3L 3 [- 

3D e-* " 3Y ^ 

X 



3L -| 
3D 9-^ 

y 



0 



(19) 



There is no definite way to discover L for any given problem, 

but it usually yields to informed guessing. 

The problem's boundary conditions present no unmanageable 

difficulties in this formulation. If 9 is specified on the 

boundary (Dirichlet boundary conditions) nothing else need 

be done. If 9 is unspecified on the boundard , the stationary 

solution will automatically fulfill the Neumann boundary 

condition, ^ =0, also known as the natural boundary condition, 
dn 

Boundary conditions of a more complicated type will require 
the functional F to be modified, but this procedure will not 
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be necessary here. It can be easily shown by using the field 
generating equations in [6], that the natural boundary condi- 
tion here results in a zero tangential electric field on the 
wire surface. 

3 . Application to the Thin-Wire Problem 

As discussed in Section II, the driving force in the 

scattering problem is the incident plane wave. Then, the 

boundary conditions to be applied are the natural condition 

along the wire, for nodes n=l, N , and a Dirichlet condition 

w 

along the outer boundary of ; h^^=h^ (total field). The 
Dirichlet condition h^^=0 has also been imposed on the top 
and bottom wire nodes . These correspond to boundary condi- 
tions along the top and bottom node rows of a larger mesh, 
as in Figure 5. 

The minimization of F for the wire scattering problem 
is carried out as follows. F is discretized by substituting 
h^j^ for which yields the approximate functional 

The exact functional F(h. ) is made stationary about the 

<po 

solution h^^ by differentiating F with respect to This 

makes F stationary about every point on Since the appro- 

ximate solution is to be found only at a finite number 
of points, the functional needs to be made stationary only at 
those points. This is done by a simultaneous solution of the 

N stationary condition equations 
w 
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0 









dh 



VRh.„ 

(f)N dh^ (})N 



dRdZ, i =N 



w 



( 21 ) 



where h. is the unknown value of h,,, at the node n=i, and 
1 (j)N ’ 



Vi^=l. Writing h^^^ explicitly, 



N c. 

/ [J VRu.'VR h u -Ru. 'y' h u ]dRdZ=0, i=l,Nw (22) 

R 1 / j n n 1 / ^ n n ’ ’ 

i,! n=l n=l 



N 

I 

n=l 



N 

or, since ^ ^ h^ is a linear operator with no coordinate 
n=l 

dependence , 






VRu -Ru.u ]dRdZ=0, i=l,Nw 
n 1 n 



(23) 



n=l 



Boundary conditions at the top and bottom wire nodes 

are already included in these equations . Application of the 

total field boundary condition is accomplished by observing 

that the h , n=N +1,N are known, where, using (15), ho=h, . 

n w — D — L 

So by splitting the sum the useful minimization equations 

for F(h^-,) are found to be 
(pN 




C. h 
in n 



+ 



N 




C. h 
in n 



n=N +1 
w 



0 



, i=l,Nw 



(24a) 
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where 



C. = / C^VRu.-VRu -Ru.u ]dRdZ (24b) 

in / R 1 n 1 n 

9. 

The matrix equation to be used for computation can be 
formed directly from (24). First, recall that the basis 
functions u are zero on and outside of the boundary 8Q . 

This means that the product of non-overlapping basis func- 
tions is zero. For example, in Figure 3, u^*u^=0. Thus 

u.u = 0 , In-il > 2 (25a) 

and 

VRu.-VRu =0 , In-il >2 (25b) 

This behavior is very useful in that it produces banded 
coefficient matrices. Special techniques can be used to 
invert large banded matrices very quickly. Writing out a 
few equations should make the matrix construction clear. 



i = l: 

i = 2 : 
i=3: 



Ciihi-HCi2h2 



^21^l'^^22^2‘^^ 



^3 2 ^ 2 ''’^ 33 ^ 3 '*’^ 



^IN +1 \ +1 ^ “^IN +2\ +2 ° 

w w w w 

23^3 ^2N +2^n +2'^‘^2N +3^ +3 

w w w w 

34^4 ^3N +3^N +3"^^3N +4^ +4 

w w w w 



0 

0 

(26) 
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The system equation can therefore be written 



Ah = -Bh„ 
w B 



(27) 



where the N xN coefficient matrix A contains the basis 
w w ~ 

function integrals for i=l,N ,n=l,N and the N xN +1 matrix 

w w w w 

B contains the integrals for i=l,N ,n=M +1,N. Both A and B 
^ w w ~ ~ 

are diagonal and have the forms 



X X 




X X 




B: 




XXX 0 




XX 0 


XXX 




X X 


0 "x X X 




0 XX 


X X 




X X 



The wire field h^ due to the boundary field hg is now 
known, since 

h =-A"^Bh„ = V h^ (28) 

— ^ ^ — D ^ — L 



Equation (28) is the forward path in the system of Figure 1. 

The actual calculation of the A. , B. will be discussed in 

in in 

Section III.C. 

B. SCATTERED FIELD INTEGRAL 

1 . Vector Potential Formulation 

The development of the feedback path is less compli- 
cated than the forward path construction given in III. A. The 
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thin-wire approximation is applied to the familiar classical 
vector potential A, and a matrix equation formulation built 
on a one-dimensional basis function expansion. The result 
will be an approximation giving the near field produced by a 
wire-surface field or current. 

From the inhomogeneous Helmholz equation 

V^A + k^A = -J (29) 



it is seen, since assumption, that A=A^Z . Then it 

follows from the vector potential defining equation 

H = V X A (30) 

that 

H, A (R,Z) (31) 

(j) 9R z 

The azimuthal invariance of H has been made explicit in (31). 

A commonly used expression for A [1] is given as an 
integration of a weighted point-source Green's function 6(r' ,r) 
over the current source where 



G(r' ,r) = 



, -jk |r-r' 
1 e o' 

4'ir T?-r ' I 



Then 



■/. 



A(r) = / J(r ' )G(r ' ,r)dV ' 

'v' 



This general situation is shown in Figure 6 



(32) 



( 33 ) 
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A(r) 




2 . Scattered Field Matrix Equation 

The formulation is easily applied to the thin-wire 
case. Recalling that h^^ = assuming that the 

current Jr, lies along the Z axis, A, for the thin wire, is 

Zj 



n^A^tR.Z) 



LI 



- j I D I 
hjz') fiTDT 



(34) 



2 2 

where, [ D | = /R +(Z-Z) , as seen in Figure 7, and the integra- 

tion is over S, the wire surface. Now, putting (34) in (31) 
yields 

R R rl/2Z 

" So 

'L 






/. 



L 

1/2 Z, 



h. (Z' )[— L ]e"2 l^ldZ' (35) 



D 
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Figure 7 



This h^^(R,Z) is approximately the scattered 

field, h^(R,Z), radiated by the wire carrying an induced 

current J„ . This wire current, in the form of h, (Z'), must 
L <po 

be approximated by a basis function expansion so that h (R,Z) 

s 

can be computed. Since h^^(Z') is a function of Z only, the 

one-dimensional linear basis functions shown in Figure 4 may 

be used. Then, for the same N nodes along the wire, as in 

w 

Section III . A. , 



N 



h, (Z') 
<po 



£ 

n=l 



h u (Z') 
n n 



(36) 



Notice that (36) defines the same nodal surface-field vector 
as (16). Substitution of (36) into (35) produces 
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hs(R,Z)=^ 



R R 




^ IDI^ |d|3 



( 37 ) 



The limits of integration reduce from ^'n-1 



Z' ,, because the basis function u (Z’) is zero for Z' ,>Z’>Z ,, 
n+1 n n-1— — n+1 



Equation (37) can be used to construct the T matrix by noting 
that (37) gives the scattered field at one point (R,Z). Then, 
by choosing a number of points in some region, a discrete 
field approximation can be constructed. The points chosen 
will be, unsurprisingly, just those Mg points along the finite 
element mesh border . So the matrix equation is 



is the coefficient of the n-th value of the wire field, and 

T, h is one term of the b-th value of the scattered field, 
bn n 

Finally, since h.j- = + — inc’ expression for the 

total field in terms of the incident and wire fields is 




(38a) 



where 





(39) 



— w 
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The scattered field h found in (38) is a near field 

— s 

expression. It is valid at any distance from the wire, but 
this behavior does not extend to the entire problem. If h^ 
were to be found at some large distance from the wire, simply 
extending the radial width of n would not suffice. This is 
because FEM calculation accuracy degrades as the distance 
from the wire to the boundary, B , increases. The magnitudes 
of the electrostatic and induction fields decrease rapidly 
as R increases, so the point is quickly reached where they 
become effectively zero as far as computer calculations are 
concerned. The resulting incorrect total field values are a 
part of the FEM R-related inaccuracies . The greater part of 
the FEM error is found in the basis function expansion of h^ . 

The planar-section basis functions used are unable to 
accurately model the near field over a large radial distance, 
where "large", in this case, may be 1/100 of the wavelength. 

To solve the far-field problem, the near-field solu- 
tion must first be obtained for the wire field h . Then (38) 

— w 

may be used to find h along any boundary. For typical far- 

s 

field problems, assumptions may be made to reduce the integral's 
complexity. Let B be the angle of the far-field point, 
measured from the Z axis at the origin. Then R=k^r sinB , 
where r is the. spherical radius of the point. The far-field 



assumption is 






|D| = k^r - Z’ cosB 


(in phase term of (38b)) 


(40a) 


n 

O 


(in denominators) 


(40b) 
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So the radiation field along a spherical arc of radius r is 



given by (38), where is 



n+1 



’bn 



sing g-jk ^ 
J 2r ^ ° 



U 



n 



(Z ' )e^^ dZ ' (41) 



n-1 



C. PROBLEM IMPLEMENTATION 

Combining equations (27) and (39) gives 

>^w = (42) 

This is, of course, the same equation as is represented in 
Figure 1, where V = -A The system "transfer matrix" 

is defined by 

h = (I - V • T)"^ • V • h. (43a) 

— w ^ ^ ^ —me 

Simple matrix algebra leads to a simplification of this 
result 

h = -(A + B • T)"^ • B • h. (43b) 

— w ^ ^ ^ —me 

Multiplying this equation by (A+BT) yields (42), so the system 

equation as derived does behave as a feedback system. 

The only variable remaining to be identified is h » 'the 

incident field vector. The incident plane wave h. can be 

^ inc 
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expressed in the normalized cylindrical coordinate system as 



h^(R,Z) = je-jZcosa 

where ' is the derivative of the zeroth-order Bessel func- 
tion. Since R will be kept small, as discussed in Section 
III.B., (44) may be simplified by replacing ' (Rsina) 
with the first few terms of its series expansion 



J ' (Rsina) 
o 



^Rsina,,_^ ^Rsina,, 3 ^Rsina^ 
^ 2 ^ ^ 16 384 ^ 



(45) 



Then rna.y be calculated at the boundary nodes by use of 

( 44 ) and ( 45 ) . 

Now that expressions for all of the quantities in (42) have 

been found, the computer program can be devised to carry out 

the calcuations . The use of an existing routine, CSMINV, for 

the inversion of (A+BT) reduces the number of major computations 

to two: equation (24b) for the A and B elements, and equation 

(38b) for the scattered field coefficients T, . 

bn 

The integration of (34b) becomes algebraicly messy, but is 
otherwide straightforward. Because of the basis function! s 
behavior, the integration area is reduced from n to , the 
region of non-zero overlap of u^ and u^. This integration 
over may then be written as a sum of integrations over the 
triangular elements comprising 3fi , since u^ is piecewise planar 
in the elements. Also, note that is independent of Z in- 
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the sense that the integration for any pair (i,n) of nodes, 
will be the same if that pair is translated along the wire. 
This means that one wire node may be selected, and the inte- 



grals concerning that node and its immediate neighbors will 

be the only ones required to construct A and B. For example, 

select i=2. Then C C 2 j^ ^2 ^2N +3 only 

w w 

integrals that need to be calculated. (Node three is not 
being neglected. Equation (24b) shows that = *^ 32 ’ 

this is just moved up one element. So ‘^23~^21*^ 

In this simple narrow mesh case, the integration of (24b) 
is required over at most three elements (^ 22 ^ 5 '^■nd the size 
and shape of all elements is the same. For many problems 
this will not be the case. For large meshes, and especially 
for skewed elements , it becomes very convenient computa- 
tionally to change from a global to a local node system. 
Instead of moving from node to node, calculating all the 
coefficients about each node one at a time, an algorithm 
is devised to move from element to element, finding the 
integrals involving the three corner nodes over each element. 
This procedure effects a considerable simplification in the 
"bookkeeping" necessary to keep track of the nodal indices 



and coordinates . 

The local node system employed here is shown in Figure 8 . 
The integrals of the combinations of local nodes (k,l) are 
found over a upper and lower element. Then the integrals 
can be placed in their proper global positions . As an example , 



35 



Upper 

Element 




Lower 

Element 



Figure 8. Local Node Indices 



consider the upper element with global nodes (2,3,N +3). The 

w 

integrations of (24b) over U are found for local node pairs 
(1,1), (1,2), (1,3), (3,1), (3, 2), and (3,3). Then the local 
node pair (1,1) integral, /(1,1) is put in A as one part in 
three of /(1, 2) is one part in two of 

w 

/(1,3)=C22 in A. Also, /(3,1)=C22» /(3,2) is part of 

W 

and /(3,3) is part of 0^2 • The remaining three local node 
integrals over U are not needed in this case, but would be 
used in larger meshes . 

For the thin mesh used here, the local node procedure 
improves only slightly the algorithm efficiency as compared 
to a global node procedure. The principle reason for its use 
was that an existing and successful subroutine was available 
to do local node integrations. This routine is named VARELA 
in the program. It requires only the local node ordering 
scheme and the coordinate (R,Z) of each local node to produce 
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a three by three matrix containing the integrations . As 
discussed above, the wire-scattering problem needs only two 
calls to VARELA to obtain the necessary values. The sub- 
routine ABLOAD then fills A and B with these integrals . 

The calculation of the coefficients T, is done by the 
familiar Gaussian quadrature method, using eight unequally 
spaced points. The interval of integration in (38b) is split 
into two segments, as in Figure 9. Then a change of variable 
is affected so that the limits of integration in each sub- 
interval are X=+l to -1. The Gaussian procedure then gives 

+ 1 n 

j f(x)dx ^ ^ ^ w^[f (x^. ) + f (-x^) ] (46) 

7-1 i=l 

where the values +x. and the weight factors w. are tabulated 
in several references such as [6]. The code used to imple- 
ment these integrations is contained in subroutine GFINT . 

Equation (42) is coded in two subroutines, CFORM and 

ZCURR. Routine CFORM produces the matrix product 

(A+BT)~^B=C, while ZCURR calculates the incident field 

vector h. according to (44) and (45), then finds the wire 

field vector h =C h. 

~w ~ ~inc 
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IV. RESULTS AND ANALYSIS 



In many other time-harmonic field problems , FEM calcula- 
tions yield good results with perhaps 20 to 30 nodes per 
wavelength. With higher node densities, say 50-75 perwavelength , 
the solution breaks down due to roundoff error in the compu- 
tations. As will be shown, this does not appear to be the 
case with the field feedback formulation. Nodal densities 
as high as 100 per wavelength have been used without diffi- 
culty. Indeed, these high densities are required in order to 
achieve good convergence of the formulation. As yet, no 
satisfactory answer for this behavior has been found. 

The only other parameter available to affect the conver- 
gence of the solution is B, , the wire-to-boundary distance. 

As mentioned in section III.B.2.., this parameter can be ex- 
pected to be critical in achieving a good field approximation. 
Computed results show this to be true. 

A great deal of effort was expended in looking for errors 
in the formulation before the solution behavior with respect 
to B^ and N^ was fully appreciated. It is still possible, no 
matter how unlikely, that some error still exists. If an 
error is present in the formulation, or code, it’s effect is 
quite peculiar. By comparison with other FEM problems, the 
possible error causes tv;o to three times the "normal" node 
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density to be needed, but otherwise seems to have no effect. 
That is, the results are satisfactory in all respects save 
node density required for convergence. 

The primary comparison case chosen was a one-half wave- 
length wire with a radius of 1/50 wavelength. This case was 
a comprom.ise forced by excessive required computing for longer 
wires and questionable integral equation comparison results 
with shorter wires. These other cases were looked at, and 
some example results are included here. Their investigation, 
however, has not been exhaustively completed due to lack of 
time . 

A. INTEGRAL EQUATION COMPARISON 

The basis used for judging the convergence of the FFM 
computation was a standard and successful type of integral 
equation formulation, Hallen's equation for scattering by 
a thin conductive cylinder. This method assumes a Z directed 
current density, and finds this v;ire current by a simultaneous 
solution of two expressions for the vector potential, A. One 
equation is just equation (33). The second is a different 
form of equation (20) 

V^A + k^A = jwe E^(r) (47) 

where E‘^(r) is the scattered electric field. Simple manipu- 
lations yield the integral equation for I(z'), C7j 
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I(Z ' )G(Z ,Z ' )dZ ' 



KZ' ) 



sink (Z-Z')dZ' 





Z(Z' ) 
Z 

o 



+ CsinkZ+D cos kZ - U(Z) (48) 

where A is an impedance per unit length, C and D are con- 
stants of integration, and U contains the incident field. 
Equation (48) is coded in subroutine INTEQ . 

Results for the 1/2 x 1/50 wavelength case are shown on 

3 

pages A-1 through A-38 of Appendix A. F calculations are 
shown as a solid line, while integral equation computations 
are displayed as the "X" curve. The parameters of each case 
are given below the graph. The magnitude plots are scaled 
in milliamperes along x, and wavelengths along y. The wire 
is to be seen lying along the y axis , and is conveniently 
demarcated by the zero end values of the integral equation 
curve. The current phase plots have the same y axis units, 
and are scaled in degrees along x. 

3 

Graphs A-1 through A-24 show F convergence with an 
increasing number of nodes. Values of N are 20,40,70 and 
100, with three plana wave incidence angles (ALPHA) used in 
each case. In these plots, the plane wave is incident from 
the upper right for ALPHA=30 and 60 degrees. ALPHA=90° is 
broadside incidence. Convergence of the current magnitude 
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• o 

is plainly seen at all incidence angles. For a=90 , the 

error is 15.5% at the wire center for N =20, but only 2.5% for 

w 

N^=70. Similarly, the phase error decreases from 18° to 1°. 

Comparing the N =70 and N =100 cases, it is seen that the 

w w 

extra nodes make very little difference in the quality of the 

solution. Convergence has been essentially accomplished at 

N^=7 0, and the additional memory and computing time for N^^=100 

are not really necessary. The convergence is less satisfying 

for incidence angles other than 90°. At ALPHA=30°, for example, 

the magnitude is not skewed quite enough, while the phase error 

increases toward the wire bottom. This situation does not 

improve as N increases . 
w 

Graphs A-25 through A-38 show FFM convergence as a function 

3 

of B, . All F calculations are m.ade with N =28. The behavior 
L w 

shown in these plots is characteristic of all cases of wire 
length and radius. For much too large, the magnitude is 
too low and the phase is off by as much as 9 0°. As Bj^ decreases, 
the magnitude rises, overshoots, and then settles back towards 
the correct, convergent value. The phase moves smoothly down 

3 

to the correct value. Notice that the F curves are nearly 
identical for DR=DZ/32 and DR=DZ/64. This implies that 
convergence has been reached by the time DR has decreased to 
DZ/32, or B =.0005A . For this case the convergence point of 
B. is actually about B =.001X , as shown in A-1 through .A-25. 

Li J_i O 

Results for several other cases of wire length and radius 
are shown in Appendix B . The one and two wavelength wire 
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calculations show that more nodes are required, assuming that 
B =.0011 is still adequate since has not changed. If 70 

Li O L 

nodes are required for a good answer at Z =X /2, then for 
two wavelengths as many as 300 nodes might be necessary. 
Similarly, well over 100 nodes could be needed for the one 
wavelength wire. The one wavelength case indicates that 

3 

this is probably true. The F solution for N =100 is percep- 

w 

tibly better than the N =70 solution, especially for off-axis 

w 

incidence angles. 

Appendix C contains graphs for two different cases: 

Aj^ = l/100A^ while or Zj^=A^/4. Graphs C-1 through C-4 

again show convergence value is no longer B. =.001X , but is 

Li O 

smaller than B. =.0005A . 

L o 

Plots C-7 through C-10 show the quarter-wavelength wire 
case. It appears from the FFM curves that convergence has 
been almost reached for N =40. This calls into question 
the accuracy of the integral equation solution. This possible 
inaccuracy should not be too surprising, however, since the 
wire is relatively thick for it's length, and the integral 
equation does much better for very thin wires. 

3 

An "error isolation" check was made on the F system by 
testing the feedback path. A fictitious sinusoidal wire 
current was generated and its radiated field calculated via 
h =T h . This field was compared to that given in reference 
[8]. The comparison is shown in plots C-11 and C-12. The 
two field calculations are equal, to the fourth significant 
digit. This indicates that the feedback path formulation and 
program coding is correct. 
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B. ERROR ANALYSIS 



Consideration of the near field behavior reveals the reason 

for the critical dependence of the solution on . The magne- 

Li 

tic field magnitude in the very near field is dominated by a 
2 

1/R term, as can be deduced from equation (37). This is 
depicted in Figure 10. Also shown are two dashed lines which 
are the linear basis function expansion of the field. It is 
clear from the figure that the approximation of the field 
deteriorates as A^^ decreases, if B^^ is held constant. Thus 
it is seen why B, =.001A is adequate for A^=.02X while 
Bj^=. 00051^ to obtain the best convergence with Aj^=.01X^. It 
is noted, however, that even at this best convergence value 

3 

of B, , the F solution is not as close to the integral equation 
L) 

solution as in the A, =.02A , B =.001X case. This fact can 

L o L o 

also be explained by consideration of the basis function 
approximation's behavior. While decreasing Bj^ improves the 
field approximation very near the wire, the approximation 
becomes worse in the neighborhood of 

It thus appears that there is a limit to the usefulness 

3 

of this formulation of the F solution to the wire-scattering 

problem. The primary difficulty lies in the inaccuracy of 

the linear basis function expansion used in the FEM calculation. 

A much better basis function may contain a term proportional 
2 

to 1/R . This would greatly reduce the sensitivity of the 
solution to Bj^. 
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Figure 10 



There is another source of error which is inherent in 
any sort of computer solution: roundoff. Simple calculations 
will be affected in the sixth or seventh significant digit by 
the use of single-precision arithmetic on the IBM-360. If 
only a few operations were required, this error could be 
ignored. But when solving large systems of equations many 
thousands of operations are required, and accuracy can be 
severely degraded. 

A rough guide to the degree of trust to be placed in a 
solution can be had in the condition number of an inverted 
matrix. This number provides an upper bound on error propa- 
gation in a system due to the inaccuracy of finite-precision 
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■^1 







arithmetic. For example, in the system X = A ^Y, a condition 
number of 1000 for A ^ would cause errors in the third signi- 
ficant digit of Y to appear as errors in the third, second or 
even first significant digit of X. 

This condition number is a measure of the similarity of 
the system equations. In the wire-scatterer case, as the 
solution points are placed closer toP-ether, the field values 
at these nodes become more indistineuishable . Hence the 
equations relating the field values at neighboring nodes 
become more nearly the same. The system solution, then, 
becomes sensitive to errors in the less significant digits. 

Condition numbers of system matrices for the cases dis- 
cussed here range from about 250 to 1000. Cases with fewer 
nodes and larger values of have lower condition numbers , 
while as N increases and/or B, decreases, the condition 
number increases rapidly. This is another reason why a good 
solution is more difficult to achieve for very thin or long 
wires. The use of a basis function more suited to the field 
behavior would allow larger values of Bj^ to be employed. This 
would lower the system condition number significantly, and 
improve the solution's accuracy and trustworthiness. Other 
than this and keeping as small as possible, there is 
nothing that can be done to avoid the problems reflected in 
the condition number. 

Even if the errors mentioned above could be eliminated, 

3 

the results of both the F and integral equation calculations 
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would still not represent the fields with complete accuracy. 

The assumptions made in both formulations are not strictly 
true for any wire of finite thickness. Current flows on 
the surface of the wire, not in the center, as assumed in 
the vector potential formulation. 

Field phase variation from side to side of the wire is not 
zero, so there is an azimuthal component of current as well as the 

/N /\ 

Z-directed current assumed. And because this b-directed 
current exists, the current is not zero at the wire ends. 

But, since essentially the same types of approximations are 
made for both solution methods , it is not possible to ascer- 
tain, from the comparisons made here, the degree to which 
either solution represents the actual fields . 
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V. CONCLUSIONS 



As a result of this investigation, the field-feedback 
formulation has been shown to be a sound technique for the 
computation of electromagnetic scattering problems. It has 

3 

been demonstrated that F calculations compare very favorably 
with results obtained by the conventional integral equation 
formulation. The success of this initial development is not 
unqualified, however. The two deficiencies previously dis- 
cussed, basis function unsuitability and excessive number of 
required nodes , could present formidable difficulties in 
problems where a priori knowledge of field behavior is not 
available, or where bodies large compared to a wavelength are 
used . 

The discrepancy between required and expected node 

density is particularly troublesome, as this fact is as yet 

unexplained. If high node densities are actually required 

by the method and are not just the result of some formulational 

error, then the attendant large matrices will limit the use- 
3 

fulness of F computations because of the great amount of com- 
puter memory required. 

In any case, the flexibility of the finite-element method 
when used in the coupled-azimuthal potential formulation 
should allow numerous heretofore difficult problems to be 
attacked. The ease with which mixed boundary conditions can 
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be met for complicated inhomogeneous scattering bodies , and 
the accuracy and simplicity of the vector potential formu- 
lation promise a wide applicability for the field feedback 
formulation . 
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APPENDIX A 
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O 

cu 
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X-SCqLE=5. OOE-04 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL = . 5: AL = 0. 02: BL = 0. 001 : RLFn = 90: NW = 20 

A-l 
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X-SCRLE=5. OOE+01 UNITS INCH. 

T-5CfiLE=5. OOE-01 UNITS INCH. 

CURRENT PHRSE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=. 5: RL=0. 02: BL=0. 001 : RLFR=90: NW=20 

A- 2 
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020 SIO 010 500 OOO 500- 010“ 510“ 020- 520“ 




X-SCRLE=5.00E-04 UNITS INCH. 

Y-5CRLE=5.00E-01 UNITS INCH. 

CURRENT MRGNITUDE :FEM/T=5QLID :INT.EQ=X 
ZL=. 5: RL=0. 02; BL=0. 001 : nLFn=60: NW=20 

A- 3 
51 




X-SCnLE=5. OOE+01 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT PHASE :FEM/T=S0LID :INT.EQ=PLUS 
ZL=.5:RL=0.02:BL=0.001: RLFA^BO : NN=20 

A-4 
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X-SCflLE=2. OOE-04 UNITS INCH. 
Y-5CRLE=5. OOE-OI UNITS INCH. 



CURRENT MAGNITUDE :FEM/T=SGLID 



INT. EQ=X 



ZL=. 5: RL = 0. 02: BL = 0. 00j^j3RLFR = 30: NW = 20 
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X-SCfiLE=5. OOE+01 UNITS INCH. 

Y-3CflLE=5.00E-01 UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=. 5: AL=0. 02: BL=0. 001 : RLFA=30 : NN= 20 

A-6 
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X-SCRLE=5. OOE-04 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ-X 
ZL=. 5: AL = 0. 02: BL = 0. 001 : RLFfl = 90 : N W = 40 

A-7 
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X-SCnLE=5. OOE +01 UNITS INCH. 

Y-3CRLE=5. OOE-01 UNITS INCH. 

CURRENT PHRSE :FEM/T=SOLID :INT.EQ=PLUS 
ZL = . 5: RL = 0. 02: BL = 0. 001 : flLFR = 90 : NW = L10 

A- 8 
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X-SCRLE=5. OOE-04 UNITS INCH. 

Y-5CRLE=5. OOE-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL = . 5: RL = 0. 02: BL = 0. 001 : RLFfl = SO: NW = 40 

A-9 
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X-SCRLE=5. OOE+01 UNITS INCH. 

T-SCRLE=5. OOE-01 UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=.5:AL = 0.02:BL = 0.001: ALFA = 60: NW = 40 

A-IO 
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-025 -020 -015 -010 -005 000 005 010 015 



z 




X-SCflLE=2.00E-04 UNITS INCH. 

Y-5CRLE=5. OOE-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SGLID :INT.EQ=X 
ZL=. 5: RL=0. 02: BL=0. 001 : flLFA=30: NW=40 

A-ll 



59 




X-SCflLE = 5. OOE^-01 UNITS INCH. 

Y-SCflLE=5. OCE-OI UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=.5;nL=0.02:BL=0.001:ALFR=30:NW=UO 

A-12 
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z 




X-SCflLE=5, OOE-04 UNITS INCH. 

T-SCRLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT. 
ZL=.5:nL=0.02:BL=0.001: RLFR =90 : NW= 70 

A- 13 
61 



040 



EO=X 




X-SCPLE=5. OOE+01 UNITS INCH. 

T-SCflLE=5. OOE-01 UNITS INCH. 

CURRENT PHASE :FEM/T=SGLID :INT.EQ=PLUS 
ZL = . 5; AL = 0. 02; BL = 0. 001 ; RLFR = 90 : N W= 70 

A-l4 
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o 




X-SCRLE=5.00E-04 UNITS INCH. 

Y-3CflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL=. 5: AL=0. 02; BL=0. 001 : ALFA =60 : NW= 70 

A-15 
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X-SCflLE = 5. OOE-^01 UNITS INCH. 

T-SCflLE=5. OOE-01 UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL = . 5: RL = 0. 02: BL = 0. 001 : AL FA = 60 : NW = 70 

A-16 
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z 




X-SCflLE=2.00E-04 UNITS INCH. 

Y-SCRLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=50LID :INT. 
ZL=. 5: AL = 0. 02: BL = 0. 001 : ALFR = 30 : NW = 70 

A-17 
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p 





X-SCflLE=5.00E+01 UNITS INCH. 

T-SCRLE=5.00E-01 UNITS INCH. 

CURRENT PHRSE :FEM/T = SOLID : I NT . EQ = PLUS 
ZL=. 5: RL=0. 02: BL=0. 001 : flLFfl=30: NW=70 

A-18 
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X-SCnLE=5.00E-04 UNITS INCH. 
Y-SCflLE=5.00E-01 UNITS INCH. 



CURRENT MRGNITUDE :FEM/T=50LID :INT.EQ=X 
ZL=0. 5: RL=0. 02: BL=0. 001 : RLFfl=90: NW=IOO 

A-19 
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X-SCflLE=5.00E+01 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT PHRSE :FEM/T=SGLID :INT.EQ=PLUS 
ZL=0. 5: RL=0. 02: BL=0. 001 : RLFR=90 : NW= 1 00 

A-20 

68 
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X-SCftLE = 5. OOE-Qii UNITS INCH. 

Y-5CflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT 
ZL = 0. 5: RL = 0. 02: BL = 0. 001 : RLFR = 60: NW- 

A-21 

69 



111 



. EQ = X 
100 




X-SCfiLE=5. OOE+01 UNITS INCH. 

Y-5CflLE=5.00E-01 UNITS INCH. 

CURRENT PHRSE : FEM/T = SQL I D : I NT . EQ^PLUS 
ZL=0. 5: RL=0. 02: BL=0. 001 : RLFR=60: NW=100 

A-22 



2H 



025 -020 -015 -010 -005 000 005 010 015 
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X-SCRLE=2.00E-04 UNITS INCH. 

Y-5CRLE=5. OOE-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T-S0LID :INT.EQ=X 
ZL = 0. 5: RL = 0. 02: BL = 0. 001 : RLFR = 30: NW=100 

A- 2 3 
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X-SCflLE=5.00E^01 UNITS INCH. 

Y-SCnLE=5.00E-01 UNITS INCH. 

CURRENT PHASE :FEM/T=S0LID :INT.EQ=PLUS 
ZL=0. 5: RL=0. 02: BL=0. 001 : nLFA=30 : NW= 1 00 

A-24 

72 




X-SCRLE=1 .OOE-03 UNITS INCH. 

T-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T-SOLID :INT.EQ=X 
ZL=.5 : RL=1/100 : DR=DZ/4 : ALFR=90 

A-25 

12 




X-SCflLE=5.00E^01 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=.5 : flL=l/100 ; DR=DZ/4 : flLFn=90 

A-26 
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X-SCflLE=5. OOE-04 UNITS INCH. 

T-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL=.5 : AL=1/100 : DR=DZ/9.8 : RLFR=90 

A-27 
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X-SCflLE=5. OOE^-01 UNITS INCH. 

Y-SCRLE=5.00E-01 UNITS INCH. 

CURRENT PHASE ;FEM/T=SOLID :INT.EQ=PLUS 
ZL=.5 : nL=l/100 : DR=DZ/9.8 : nLFfl=90 

A-28 
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X-SCOLE=5.00E-0ii UNITS INCH. 



T-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL = .5 : AL=l/50 : DR^=^gZ/2 : ALFfi = 90 
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X-SCRLE=5.00E+01 UNITS INCH. 

Y-SCFILE = 5. OOE-01 UNITS INCH. 

CURRENT PHRSE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=.5 : RL=l/50 : DR=DZ/2 ; RLFR=90 

A-30 
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X-SCflLE=5.00E-04 UNITS INCH, 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SGLID :INT.EQ=X 
ZL=.5 : AL=l/50 : DR=DZ/4 : flLFA=90 

A-31 
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X-SCflLE=5.00E+01 UNITS INCH. 
Y-SCRLE=5.00E-01 UNITS INCH. 



CURRENT PHRSE :FEM/T=SOLID 
ZL=.5 : RL=l/50 : DR=DZ/4 

A-32 



INT. EQ=PLUS 
RLFn=90 
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X-SCflLE=5. OOE-04 UNITS INCH. 

T-SCRLE=5.00E-01 UNITS INCH. 

CURRENT MRGNITUDE :FEM/T=S0LID :INT.EQ=X 
ZL=.5 : nL=l/50 : DR=DZ/8 : RLFfl=90 

A-33 



ai 




X-SCRLE=5.00E+01 UNITS INCH. 

Y-SCRLE=5.00E-01 UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ-PLUS 
ZL=.5 : RL=l/50 : DR=DZ/8 : RLFA=90 

A-34 
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X-SCflLE=5.00E-04 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MfiGNITUDE :FEM/T=SGLID :INT.EQ=X 
ZL=.5 : RL=l/50 ; DR=DZ/32 : RLFR=90 

A-35 
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X-SCRLE = 5.00E^-01 UNITS INCH. 

Y-5CRLE=5. OOE-01 UNITS INCH. 

CURRENT PHASE : FEM/T = SGL I D •: I NT . EQ = PLUS 
ZL=.5 : flL=l/50 : DR=DZ/32 : flLFR=90 

A-36 
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X-SCRLE=5. OOE-04 UNITS INCH. 

Y-SCflLE=5.00E-01 UNITS INCH. 

CURRENT MfiGNITUDE :FEM/T=SGLID ;INT.EQ=X 
ZL=.5 : RL=l/50 : DR=DZ/64 : flLFfl=90 

A-37 
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X-SCflLE=5.00E+01 UNITS INCH. 

T-SCflLE=5.00E-01 UNITS INCH. 

CURRENT PHASE :FEM/T=S0LID :INT.EQ=PLUS 
ZL=.5 : RL=l/50 : DR=DZ/6^ : ALFR=90 

A- 3 8 
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APPENDIX B 



Z 



X 







X-SCRLE=5. OOE-04 UNITS INCH. 
Y-SCRLE=2.00E^-00 UNITS INCH. 



CURRENT MRGNITUDE :FEM/T=SGLID :INT.EQ=X 
ZL=2. 0: RL=0. 02: BL=0. 001 : fiLFR=90: NN=70 
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X-SCflLE=5.00E^OI UNITS INCH. 

T-SCRLE=2.00E^O0 UNITS INCH. 

CURRENT PHRSE :FEM/T-SOLID :INT,EQ=PLUS 
ZL=2. 0: RL=0. 02: BL=0. 001 : RLFR= 90 : NN= 70 
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X-SCfiUE = S.OOE-Ol| UNITS INCH. 
T-SCflUE=2.00E»00 UNITS INCH. 



CURRENT MAGNITUDE :FEM/T=SaLIO :INT.EQ=X 
ZL=2. 0: RL=0. 02: BL=0. 001 : ALFR=S0 : NN=70 



o n 




X-SCflLE=5. OOE+01 UNITS INCH. 
T-SCnLE=2, OOE+00 UNITS INCH. 



CURRENT PHRSE :FEM/T=SOLID :INT. 
ZL=2.0:flL=0.02:BL=0.001:RLFP=60: 



EQ=PLUS 
NW = 70 
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X-SCflLE=5.00E-04 UNITS INCH. 
Y-SCRLE=2.00E+00 UNITS INCH. 



CURRENT MRGNITUDE :FEM/T=SQLID :INT.EQ=X 
ZL = 2. 0: RL-O. 02: BL = 0. 001 : RLFR = 30: NW = 70 
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I 




X-SCflLE=5.00E+01 UNITS INCH. 

Y-SCflLE=2.00E+00 UNITS INCH. 

CURRENT PHRSE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=2, 0: RL=0. 02: BL=0. 00 1 : RLFR=30 : NW=70 
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X-SCRLE=5. OOE-OU UNITS INCH. 
T-SCRLE= I . OOE^OO UNITS INCH. 



CURRENT MAGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL=1 . 0: RL=0. 02: BL=0. 001 : RLFfi=90: NW=70 
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X-SCflLE=5. OOE+01 UNITS INCH. 

T-SCflLE=l .OOE+00 UNITS INCH. 

CURRENT PHRSE : FEM/T = SOL I D •: I NT . EQ-PLUS 
ZL=1 . 0: RL=0. 02: BL=0. 001 : RLFR=90: NW=70 
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X-SCRLE=5. OOE-04 UNITS INCH. 
Y-SCflLE=l.OOE^CO UNITS INCH. 



CURRENT MRGNITUDE :FEM/T = SGLID :INT.EQ=:X 
ZL=1 . 0: flL=Q. 02: BL=0. 001 ; RLFR=60 : NN=70 
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X-SCRLE=5.00E+OI UNITS INCH. 

Y-3CRLE=l .OOE^-00 UNITS INCH. 

CURRENT PHASE :FEM/T=SQLID :INT.EQ=PLUS 
ZL=1 . 0: RL=0. 02: BL=0. 001 : RLFR=60: NN=70 
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X-SCRLE=5,00E-04 UNITS INCH. 
T-3CnLE=l .OOE+00 UNITS INCH. 



CURRENT MRGNITUDE :FEM/T=SOLID :!NT.EQ=X 
ZL = 1 . 0: RL = 0. 02: BL = 0. 001 : flLFR = 30 : NW = 70 
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X-SCRLE=5.00E^01 UNITS INCH. 
'r-SCRLE=l . OOE + 00 UNITS INCH. 



CURRENT PHASE :FEM/T=SOLID :1NT.EQ=PLUS 
ZL=1 . 0: AL = 0. 02: BL = 0. 00 1 : RLFR = 30: NW = 70 



98 



riOO coo 200 TOO OOO TOO- 200- coo- hOO- 500- 



I 




X-SCflLE=5. OOE-04 UNITS INCH. 
T-SCRLE=I.00E+00 UNITS INCH. 



CURRENT MAGNITUDE :FEM/T=SOLID 



ZL=1 . 0: RL=0. 02: BL=0. 001 : RLrR=90 



: INT. EQ = X 
: NN= 100 
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X-SCflLE=5.00E+0I UNITS INCH. 
T-5CflLE=l.00E^-00 UNITS INCH. 



CURRENT PHRSE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=1 . 0: RL=0. 02: BL=0. 001 : RLFR=90 : NW= 1 00 
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I 




X-SCRLE=5.00E-04 UNITS INCH. 

T-SCRLE=l . OOE-^00 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID :INT,EQ=X 
ZL=1 . 0: AL=0. 02: BL=0. 001 : RLFn=60: NW=IOO 
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X-SCRLE=5. OOE+OI UNITS INCH. 
Y-SCRLE=1 .OOE^OO UNITS INCH. 



CURRENT PHRSE : FEM/T = SOL I D • : I NT . 
ZL=1 . 0: RL=0. 02: BL=0. 001 : RLFR=60: 



EQ=PLUS 

NW=100 



102 




X-SCRLE = 5.00E-0Li UNITS INCH. 

T-5CRLE=l .OOE^CO UNITS INCH. 

CURRENT MRGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL=1 . 0: RL=0. 02: BL=0. 001 : flLFR=30: NN=IOO 
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X-SCflLE=5.00E+01 UNITS INCH. 

T-SCflLE=l . OOE+QO UNITS INCH. 

CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=1 . 0; AL = 0. 02: BL = 0. 001 : RLFA = 30 : NN= 1 00 
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X-SCflLE=5.00E-04 UNITS INCH. 
Y-SCRLE=5. OOE-01 UNITS INCH. 



CURRENT MAGNITUDE : FEM/T=SOL I D :INT.EQ=X 
ZL=. 5: AL=0. 01 : BL=0. 001 : RLFR=90 : NW=40 
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X-SCflLE=5. OOE+01 UNITS INCH. 
Y-SCflLE=5.00E-01 UNITS INCH. 



CURRENT PHRSE :FEM/T=SOLID :INT.EQ=PLUS 
ZL = . 5: RL = 0. 01 : BL = 0. 001 : flLFfl = 90 : NW = 40 
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X-SCRLE=5.00E-04 UNITS INCH. 
Y-5CflLE=5.00E-01 UNITS INCH. 



CURRENT MAGNITUDE :FEM/T=S0LID :INT.EQ=X 
ZL=.5:AL=0.01: BL=0 . 0005 : RLE A=90 : NW=40 
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X-SCRLE=5.00E^01 UNITS INCH. 
T-SCflLE=5.00E-01 UNITS INCH. 



CURRENT PHASE :FEM/T=SQLID :INT.EQ=PLUS 
ZL=. 5: RL = 0. 01 : BL = 0 . 0005 : nLPn = 90 : NN = 40 
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X-SCflLE=5.00E-04 UNITS INCH. 
Y-SCflLE=5. OOE-01 UNITS INCH. 



CURRENT MAGNITUDE : FEM/T==SOL I D :INT.EQ = X 
ZL=. 5: RL=0. 01 : BL=0 . 0005 : ALFR=90 : NW=70 
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X-SCflLE=5. OOE+01 UNITS INCH. 
Y-SCRLE=5.00E-01 UNITS INCH. 



CURRENT PHASE :FEM/T=S0LID :INT.EQ=PLUS 
ZL=. 5: RL=0. 01 : BL=0 . 0005 : fiLFA=90 : NN= 70 
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X-SCRLE=5. OOE-OU UNITS INCH. 

T-5CnLE=5.00E-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID ;INT.EQ=X 
ZL=. 5: AL=0. 01 : BL=0 . 0005 : RLFR=90 : NW= 1 00 
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I 




X-SCRLE=5. OOE+01 UNITS INCH. 
Y-SCRLE=5.00E-01 UNITS INCH. 



CURRENT PHASE : FEM/T=SOL I D :INT.EQ=PLUS 
ZL=. 5: AL=0. 01 : BL=0 . 0005 : nLFR=90 : NW=100 



112 



800 000 fiOO ^0O 000 200- hOO- OOC- 900- 010' 



z 




X 



X-SCflLE=l . OOE-04 UNITS INCH. 

T-SCRLE=2. OOE-01 UNITS INCH. 

CURRENT MAGNITUDE :FEM/T=SOLID ;INT.EO=X 
ZL=. 25: RL=0. 01 : BL=0. 001 : ALF R= 90 : NW =40 
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X-SCflLE=5. OOE+01 UNITS INCH. 
T-SCnLE=2.00E-01 UNITS INCH. 



CURRENT PHASE :FEM/T=SOLID :INT.EQ=PLUS 
ZL=. 25: RL = 0. 01 : BL = 0. 001 : RLFfl = 90 : NW = 40 
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X-SCflLE=1.00E-04 UNITS INCH. 

Y-SCflLE=2.00E-01 UNITS INCH. 

CURRENT MRGNITUDE :FEM/T=SOLID :INT.EQ=X 
ZL=. 25: nL=0. 01 : BL=0. 0005; RLFfl=90: NW-70 
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X-SCnLE=5, OOE+01 UNITS INCH. 
T-SCnLE=a. OOE-01 UNITS INCH. 



CURRENT PHRSE :FEM/T=SOLID :INT.EQ-PLUS 
ZL=. 25; RL=0. 01 : BL=0. 0005: RLFR=90 : NW-70 
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X-SCflLE=5. OOE+02 UNITS INCH. 
T-SCflLE=5.00E-01 UNITS INCH. 



BOUNDARY FIELD MAG. : ZL= 
SOLID^EXACT : DIAMONDS=T 



.5 : DR=DZ/4 
: AL=1/100 
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X-SCflLE=l .OOE+00 UNITS INCH. 

T-SCRLE=5.00E-01 UNITS INCH. 

BOUNDRRT FIELD PHRSE : ZL=.5 : DR=DZ/4 
SOLID=EXRCT : TRIRNGLES=T : RL=1/100 
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